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I. INTRODUCTION 

The finite size analogs of two-dimensional (2D) Wigner 
crystals have received considerable attention from theo- 
reticians 0, aSiiilliili3IIl|lllil|ll 
llSlllGl. and more recently also from experimental groups 
Particles in mutual repulsion can be 
confined using a variety of methods, most importantly 
electrostatic traps of the Paul or Penning types, and sim- 
ple hard walls, usually of circular shape. The repulsion 
between particles depends upon the details of the ex- 
periment. Arrays of electrons on the surface of liquid 
helium or in quantum dots [2^ can be described 
by conventional Coulomb forces. However, logarithmic 
forces are more relevant for vortex clusters in superfluids 
[2^ or Bose-Einstein condensates [S^. Mesoscopic col- 
loidal particles undergo dipolar interactions within 
an external magnetic field j26| |. or a Yukawa (screened 
Coulomb) potential when they are charged and placed in 
a viscous medium. 

Finite 2D assemblies have also been synthesised on a 
more macroscopic scale, using metal particles close to 
1mm in size (19]. In this study, the authors were able 
to visualise and locate not only stable configurations, 
but they could also estimate the shape of the transition 
states, which are here defined as stationary points with 
Hessian index one |4^ . 

Another example of experimentally accessible Wigner 
crystals occurs for dusty plasmas. These highly charged 
particles of micrometer size and charges up to Z = 10^ 
can be stabilised as gaseous plasmas. The gravitational 
forces and the electric field balance each other, so that the 
dust particles form highly regular two-dimensional struc- 
tures. By changing the electric field, or conversely the 
pressure, of the plasma, transitions from solid to liquid- 



like phases are observed A two-stage melting was 

proposed, which does not follow the predictions of two- 
dimensional melting theories, but rather goes through an 
intermediate phase. This phase has been described as an 
oscillating crystal . Finally, by measuring the spectra 
of dust particles, it has been suggested that the dynamics 
may be a size-dependent phenomenon, that is the effects 
of various types of motion such as intershell rotation and 
breathing may play different roles |2^ . 

A major part of our theoretical understanding of the 
above systems comes from numerical studies. The static 
properties of the assemblies were obtained using opti- 
mization methods ranging from simulated annealing to 
genetic algorithms, and a Mendeleev- type table was pro- 
posed by several authors Ij'F] to reflect the special stabil- 
ity of certain sizes. Comparisons with the theoretical pre- 
dictions of the 'Thomson atom' [23 were also carried out 
0. Shell effects were more generally investigated semi- 
analytically by Koulakov and Shklovskii • Metastable 
configurations (excited states) and the potential energy 
barriers between them have been studied in Refs. 0, |l2j . 
[T5I I and [S^l- Topological defects were also identified for 
stable structures (l3j . 

Melting in 2D Coulomb clusters has been a subject of 
debate. The results of the classical and quantum Monte 
Carlo simulations of Lozovik and Mandelshtam (^ were 
interpreted as a two-step melting process, where orien- 
tational melting occurs prior to radial melting. In this 
process, shells remain quasi-crystalline, but still undergo 
some concerted rotational motion. Radial melting con- 
versely involves exchanges between particles of adjacent 
shells. This two-step process is generally accepted for 
small assemblies, but it is not clear whether, in large 
clusters, the external shell may also exhibit orientational 
melting |T5| or not 0. 
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Several factors are likely to influence the stability of 2D 
Coulomb clusters. First, rotational barriers are strongly 
dependent on the structure, which in turn depends on the 
size. More stable structures are obtained when the core 
is based on a triangular lattice, thus resembling most 
observed Wigner crystals. More importantly, the very 
different nature of orientational and radial melting has 
led to some confusion about how suitable observables 
should be defined to probe these processes separately 
[T^ l2^ [29ll. In their study of binary 2D clusters, Drocco 
et al. plq chose an order parameter based on the polar 
angles of the particles in order to quantify the extent 
of orientational motion in Langevin molecular dynamics 
simulations. However, this measure becomes impractical 
as soon as a central particle is present, as it then becomes 
very noisy. 

Melting in atomic and molecular clusters also dis- 
plays a very rich phenomenology |f30l |. In particular, sur- 
face melting jsj, plastic transitions js^l or anomalously 
high melting temperatures [s^l have been identified 
from experiments or simulations. For three-dimensional 
trapped-ion clusters, the structures are sometimes re- 
lated to those found in atomic clusters [s^] ■ The study of 
atomic clusters also shares with Coulomb assemblies the 
lack of a proper universal order parameter to characterise 
the melting point. For instance, the Lindemann criterion 
for melting does not always provide useful information in 
the case of solid-solid transitions and premelting effects 

M 

The largest Lyapunov exponent A, which quantifies 
chaos and the sensitivity to initial conditions, is a use- 
ful tool for probin g th e complex dynamics of atomic and 
molecular clusters |3a 123, lla |33, , especially at melt- 
ing. Because Lyapunov exponents are related to the sec- 
ond derivatives of the potential energy surface (PES) [4ll| . 
we expect them to be sensitive to orientational melting. 

In this paper, we revisit the dynamics and thermody- 
namics of 2D Coulomb clusters by focusing on dynamical 
measures such as A or the diffusion coefficient, D, and by 
performing extensive Monte Carlo simulations employing 
the parallel tempering strategy 0| . We also consider the 
interplay between the dynamics and the energy landscape 
of these clusters by constructing disconnectivity graphs 
[4^ 13 ■ Finally, we investigate how the thermodynamic 
behaviour evolves with size. 

This article is organised as follows. After briefly pre- 
senting our methods, we choose some specific examples 
to illustrate the difficulties associated with an unambigu- 
ous definition of the melting point in these clusters, and 
we then discuss the large size variations. We finally sum- 
marise and conclude in Section IV. 



II. METHODS 

2D Wigner clusters are characterised by the following 
classical potential energy: 

^(R) = E fl^ exp(-^r,,) -f Arl (1) 

In this equation, {qi} = 1 or 2 are the charges carried 
by the particles, Tij is the distance between particles i 
and j , and is the distance from particle i to the clus- 
ter centre-of-mass. The parameter p will usually be 1 
(Coulomb case), but we also specifically studied the case 
p = 3 (dipolar interaction) for which some results will 
be selected below. The range of the electrostatic po- 
tential, K, was set to (pure Coulomb) or to a finite 
value (shielded Coulomb). Finally,, and for comparison 
with the work reported in Ref. on binary clusters, 
we have taken A = 10. It should be noted that, in 
the absence of the exponential term, the physics of the 
system is independent of A: under the transformation 
(r, E) (rAi/(2+P)^ £;/^P/(2+p))^ the energy can be cast 
in a form that does not depend on A. In the following 
presentation, unit masses and charges are used and all 
properties are explicitly given in dimensionless reduced 
units. 

The global minima of the clusters were investigated 
using the Monte Carlo plus minimization algorithm, also 
known as basin-hopping . The 2D clusters can 

generally be described as multishell systems, and we will 
use the (ni,n2...) notation for a homogeneous cluster 
having Ui particles in its shell. Binary clusters will 
be referred to as {Ns^Nd}^ where and Nd are the 
respective numbers of singly- and doubly-charged parti- 
cles. Due to the relatively long range of the potential 
global optimisation is relatively straightforward for 
Coulomb systems. 

Beyond the stable isomers, we also investigated higher- 
index stationary points on the energy landscape. Locat- 
ing the saddle points with Hessian index one enabled us 
to construct disconnectivity graphs '47^ , in order to relate 
the observed dynamical and thermodynamical behaviour 
to the potential energy surface for these clusters. We 
refer the reader to Ref. for further details of how sta- 
tionary points are located and disconnectivity graphs are 
generated. 

The dynamics has also been studied by solving Hamil- 
ton's equations of motion at constant energy. The cou- 
pled equations were integrated by the 4-step Runge- 
Kutta method, which keeps the energy constant to at 
least seven digits. For the first trajectory, the initial 
positions of the particles were obtained from the global 
minimum structure and a random set of momenta were 
chosen. The momenta were then scaled to set the lin- 
ear and angular momentum to zero and the vibrational 
temperature to around 10~^. The time step was taken 
to be 0.0025. The first 10^ steps were used to equili- 
brate the system and discarded, then 10^ steps were used 
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for each simulation. Further trajectories start with the 
configuration obtained from the previous run, with ve- 
locities scaled approximately to give a new temperature 
increased by a factor 1.25. Again, the equilibration steps 
were discarded. 

For each trajectory the average kinetic temperature, 
root-mean-square (RMS) bond length fluctuation index 
(5, the maximal Lyapunov exponent. A, (MLE) and diffu- 
sion coefficient, D, were calculated. 6 is commonly used 
to detect melting in bulk or finite (3D) atomic systems, 
as it is expected to jump from a low value to around 0.1- 
0.15. Both the diffusion coefficient and the Lyapunov 
exponent were calculated from shorter parts of the entire 
trajectory. D was obtained from the slope of the mean 
square atomic displacement. This measure describes the 
breathing and inter-shell exchanges (isomerizations) , but 
does not identify even qualitative changes in the rota- 
tional structure. For homogeneous systems it is very dif- 
ficult to define measures of the angular motion. Since 
the total angular momentum is zero, one can only mea- 
sure angular velocity and momenta for separate shells; 
however, the shell definitions lose their significance once 
isomerizations start taking place. In addition, the an- 
gular properties become ill-defined for clusters having a 
central particle. 

The calculation of the MLE was motivated by sev- 
eral studies indicating possible relations between phase 
changes and this dynamical indicator. The MLE mea- 
sures the sensitive dependence on initial conditions and 
is defined from the exponential separation of two tra- 
jectories that start infinitesimally close together in the 
phase space. A tangent space approach was used in the 
numerical calculations |49 | . However, we should mention 
that semi-analytical theories have also been proposed to 
quantify the MLE from the statistical properties of the 
potential energy surface [49|. 

The equilibrium thermodynamics of the 2D clusters 
was also studied by means of parallel tempering Monte 
Carlo (MC) simulations. Parallel tempering is a par- 
ticularly efficient method to accelerate convergence in 
systems whose natural dynamics is affected by broken 
ergodicity. Most often such problems are related to mul- 
tiple funnels or 'glassy' energy landscapes. The MC 
simulations were performed in the canonical ensemble, 
with temperatures in the range 10~^ < T < 10. About 
80 runs were carried out in this range, with the tem- 
peratures nearly equally spaced on a logarithmic scale. 
10^ MC cycles were used for equilibration, followed by 
10^ MC cycles for the actual calculations. The MC cal- 
culations were principally used to compute the thermal 
properties, such as the heat capacity, Cy 



III. RESULTS 

The snapshots obtained from experiments on colloidal 
particles [13 or dust plasmas [23, |21| , as well as the vari- 
ous molecular dynamics (MD) simulations, exhibit three 
types of characteristic motion. At low energy or tem- 
perature, or for dust plasmas under high pressures, only 
small amplitude motion is observed, in which the parti- 
cles oscillate around their mean positions. As the energy 
increases, these oscillations gain enough momentum to 
achieve complete rotations in their respective shells. The 
directions of these rotations are random, and they may 
change in time. However, in order to preserve the total 
angular momentum, different shells must rotate in op- 
posite directions. Eventually, actual isomerizations be- 
tween shells may take place, as previously examined in 
3D systems [s^- In one case, it has also been reported 
that isomerizations may occur before the full rotational 
motion is achieved We expect this observation to 

be due to rather infrequent sampling of the experimen- 
tal trajectories ^2^. It is clear that the angular and the 
radial motion of particles make different contributions to 
all the measures aforementioned. 




FIG. 1: Dynamic and thermodynamic properties of some ho- 
mogeneous 2D Coulomb clusters with n = 17, 21, 22, and 
30 as a function of temperature, (a) Maximum Lyapunov 
exponent (MLE) A; (b) RMS bond length fluctuation S; (c) 
diff'usion coefficient D; (d) canonical heat capacity Cv The 
data in (a) to (c) are from microcanonical MD simulations, 
the data in (d) are from canonical Monte Carlo simulations. 
C'v is in units of the Boltzmann constant 



In both MC and MD simulations the relation between 
the energy and the temperature is quasi linear. 



In order to study the mechanism of phase transitions 
in 2D clusters, we have carried out simulations of homo- 
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geneous clusters systems in the size range 3 < n < 40. 
We selected several sizes displaying different but typical 
characteristics. In Fig.^ ^, D, S and Cy are plotted with 
respect to the average temperature for sizes 17, 21, 22 and 
30. The global minima for these clusters are found to 
be (1,6,10), (1,7,13), (2,842) and (5,10,15), respectively, 
in agreement with Ref. 0- The MLE starts becoming 
nonzero at very low temperatures, around T = 10~^. In 
this regime, the dynamics are more irregular for n = 17 
and 21. Similarly, the Lindemann index and the diffu- 
sion coefficients are also higher for n = 17 and 21 in the 
low T regime. Both these clusters exhibit an interme- 
diate phase over a broad temperature range (note the 
logarithmic temperature scale). On the other hand, the 
30-particle cluster exhibits a single phase transition. 

While the three dynamical indicators agree with one 
another, the heat capacities provide rather little infor- 
mation in comparison. The general behaviour observed 
for all sizes is that of harmonic oscillators at very low 
temperatures {Cy — > (2n — 3)kB) or for free, indepen- 
dent particles at high temperatures (Cy — > (2n — S)kB 
again due to confinement). In the intermediate temper- 
ature range, a strong decrease is seen near T ^ 2 for all 
sizes, following an order-disorder peak below T = 0.1. 
The curves for a specific cluster are built on the basis of 
this generic behaviour, which reminds us of the caloric 
curves computed for Lennard- Jones polymers |50l |. For 
some sizes, such as n = 21 or 30, extra bumps or shoul- 
ders occur as a consequence of isomerizations. In general, 
the strong increase of A, D or (5 at T ~ 0.03 correlates 
with a peak in the specific heat. However, the phenom- 
ena responsible for the peculiar dynamical features for 
the 17- and 21-particle clusters do not have any thermal 
signature. The premelting heat capacity peak seen for 
n = 21 at T = 0.01 does not have any counterpart in the 
dynamical indices either. 

These data can be interpreted using the disconnectiv- 
ity graphs shown in Fig. |21 In this figure, saddles con- 
necting two permutations of the same isomers have been 
included to account for the rotational barriers. Such sad- 
dles are manifested by the presence of 'copies' of the cor- 
responding minima on the graph. 

The energy barriers, which link the global minimum 
to itself through an internal rotation, are much lower for 
n = 17 {AE = 4 X 10"*^) and n = 21 {AE = 10"^, this 
is nearly a case of free rotation) than for larger clusters, 
where they are closer to AE = 10~^. Such degenerate 
rearrangements ^3 do not play a major role in thermal 
equilibrium, but can have important dynamical conse- 
quences. This is precisely what we observe in Fig. ^ 
The barriers for rearrangements are much higher in en- 
ergy than those involved in the rotational motion, but in 
some cases other stable isomers lie relatively close to the 
global minimum. For n = 21 and n = 30, such isomer- 
izations are indeed seen and correlate with bumps in the 
canonical heat capacity. 

^From the previous results it is clear that one cannot 
draw general conclusions by looking at dynamical indi- 



cators alone, or at observables aimed at characterising 
thermal equilibrium. We have also looked at two specific 
cluster sizes, for which the results of the pure Coulomb 
case have been compared to other interaction forms. The 
31-particle cluster was studied with both 1/r and in- 
teractions, the latter expression being more relevant to 
colloidal particles, such as the ones involved in the ex- 
perimental setup of Ref. il^ The 42-particle cluster has 
recently been the focus of an experimental report on dust 
plasmas In this situation, a screened Coulomb or 

Yukawa form e~"' /r is more appropriate than the bare 
Coulomb law. Here we chose k = 1, but the results pre- 
sented hereafter in Fig. O are not significantly different 
for other values in the range 0.1 < k < 10. 

We first discuss the 31-particle clusters, which display 
different behaviour depending on the form of the inter- 
action. While the Coulomb potential leads to a (5,11,15) 
3-shell global minimum, the shorter range dipolar inter- 
action favours a (1,5,10,15) ground state. Because this 
structure has a central particle, the vibrational motion is 
more constrained, and the harmonic entropy of (5,11,15) 
is actually larger. This competition between energy and 
entropy leads to a structural transition, which is seen 
in the heat capacity of Fig. |3fa) as a significant peak 
near T — 0.03, while the melting peak is located at 
T 0.5. The Coulomb cluster behaves very differently, 
with a single Cy peak at T = 0.03. Scrutiny of the Lin- 
demann index indicates that this cluster undergoes free 
internal rotation at very low temperatures T < 10"^, 
before fully melting at T = 0.03. In the case of the dipo- 
lar interaction, only a single, but very sharp, increase is 
seen, at temperatures even higher than the melting peak, 
T ^ 0.1. This surprising result has been checked and con- 
firmed by extensive sampling. After regular quenching of 
the MD trajectories in this temperature range, we found 
that the sharp increase of d was correlated with a sudden 
isomerization to the (5,11,15) second-lowest minimum, 
which is also involved in the preliminary structural tran- 
sition. The lowest energy transition state from the global 
minimum was found to be significantly higher than the 
energy of the second isomer. In such a situation, MD sim- 
ulations suffer from broken ergodicity, since they cannot 
sample the available phase space until sufficient energy 
is added. In cluster physics, a well-studied example of 
broken ergodicity is provided by the 38-atom Lennard- 
Jones cluster |5ll [S^- The present sharp jump of the 
Lindemann index should not be interpreted as the onset 
of melting, but rather as a signature of isomerization. 

42-particle clusters were investigated with the aim of 
finding some thermodynamic confirmation that melting 
proceeds in multiple steps. Unfortunately, neither the 
specific heats of Fig. IHIb) nor the Lindemann indices of 
Fig. |3{d) show any multiple-stage behaviour that could 
be related to the experimental observation ^2l'| . We thus 
conclude that the dust plasma system studied by Melzer 
and coworkers cannot be simply modelled as if it were 
made of a single layer. 

Comparing the general curves for the Coulomb, dipolar 
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FIG. 2: Disconnectivity graphs for the homogeneous clusters of Fig. The global minimum energy is shifted to zero for each 



and Yukawa potentials also provides useful information 
about the factors influencing the melting point. Yukawa 
systems exhibit a significant shift of their melting tem- 
perature toward lower T . Since the Yukawa potential is 
an attenuated form of the pure Coulomb form, this is the 
expected trend. At large distances, the dipolar interac- 
tion is also much smaller than the Coulomb interaction. 
However, the melting temperature of Coulomb clusters 
is about one order of magnitude smaller than for the 
dipolar cluster. Hence the long-range part of the poten- 
tial plays opposite roles in the dipolar and Yukawa forms. 
This result suggests that the location of the melting peak 
is mainly driven by the short-range repulsive part of the 
potential. 

We now turn to small, binary Coulomb clusters, such 
as those recently investigated by Drocco and coworkers 
[T^ . Ns particles with charge 1 and Nd particles with 
charge 2 are stabilised in a parabolic trap with A — \Q. 
The structures again consist of concentric shells where 
the singly charged particles form the inner shells in or- 
der to minimize repulsion. Drocco et al. used Langevin 
MD simulations for mixed clusters containing Ng = 3-34 
singly charged and Nd = 5,6 and 7 doubly charged parti- 
cles, and they also determined the melting points. These 
authors concluded that the melting behaviour of these bi- 
nary trapped particles strongly depends on whether the 
numbers of particles in the inner and outer shells are 
commensurate. They also found that the highest melt- 
ing points were obtained for Ng — Nd + 1, where the 
corresponding structures all have a single particle at the 
centre. 

We have performed MD simulations for these hetero- 
geneous systems under similar conditions to those of the 
homogeneous clusters. In Fig. ^ A, Z?, 5 and the spe- 



cific heat Ci, are plotted for {4,5}, {5,5}, {6,5} and 
{7,5}, selected as further typical examples, highhght- 
ing differences in the stability and dynamics. The on- 
set of chaos occurs at temperatures close to those seen 
in the homogeneous systems. {6,5} is the most sta- 
ble cluster, in the sense that it has the highest melt- 
ing point in Ref. 0, while {5,5} is the most irregular. 
Both the diffusion and Lindemann indices show the same 
qualitative trends, namely that the melting points, as 
measured by the strong increase in (5, follow the order 
{4,5}<{7,5}<{5,5}<{6,5}. 

As for homogeneous clusters, the variation of the heat 
capacities does not exactly match the dynamical indices. 
Here the generic shape again starts at a constant value 
at low temperature, which corresponds to harmonic be- 
haviour. The high temperature gas phase is still char- 
acterised by the same constant value. In the interme- 
diate temperature range, two main peaks are associated 
with the melting of successive subclusters made of singly 
charged, then doubly charged particles. Nearly one order 
of magnitude separates the two associated melting tem- 
peratures. This interpretation was further confirmed by 
calculating the heat capacity of the ternary {7,6,6} clus- 
ter made of 7 singly charged particles, 6 doubly and 6 
triply charged particles. For this cluster, three clear fea- 
tures are observed in Ct,, and only above T ~ 2 do the 
triply charged particles mix with the rest of the cluster. 

On top of this generic behaviour, a specific isomeriza- 
tion shoulder is seen for the heat capacity of the {7,5} 
cluster in Fig. ^d), at the same temperature T ~ 0.01 
where the three dynamical parameters increase signifi- 
cantly. Again, the free rotation of some internal shells 
revealed by A, D and 5 at low temperature has no ther- 
mal signature. The disconnectivity graphs for these clus- 
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FIG. 3: (a) and (b): heat capacities of 31- and 42-particle 
clusters, respectively (in units of the Boltzmann constant); (c) 
and (d): RMS bond length fluctuation indices. For n — 31, 
the results of the Coulomb (1/r) and dipolar (l/r^) interac- 
tions are shown. For n = 42, the results obtained with the 
Yukawa potential {e~^ /r) are compared to those with the bare 
Coulomb form. 



FIG. 4: Dynamic and thermodynamic properties of some bi- 
nary {n,5} 2D Coulomb clusters with n = 4-7 versus tem- 
perature, (a) Maximum Lyapunov exponent (MLE) A; (b) 
RMS bond length fluctuation S; (c) diffusion coefficient D; 
(d) canonical heat capacity Cv The data in (a) to (c) are 
from microcanonical MD simulations, the data in (d) are from 
canonical Monte Carlo simulations and Cv is in the units of 
Boltzmann constant. 



ters are shown in Fig. [S] As is particularly obvious from 
the graphs for {5,5} and {6,5}, commensurate structure 
between the singly- and doubly-charged shells is not the 
only factor needed to explain the complicated behaviour 
exhibited by these finite Wigner crystals. The extremely 
low barriers for the {4,5} and {5,5} clusters permit facile 
rotation, which is correlated with the rapid increase of 
all dynamical indicators at low temperatures. The {6,5} 
cluster is especially resistant to isomerization because an 
extra, central particle hinders the radial motion of the 
particles in the other shells. Only for the {7,5} cluster 
is the isomerization energy low enough to give rise to 
an extra peak in the heat capacity. This case, and the 
case of the single-isomer {6,5} cluster, further show that 
the multiple-peak structure of the heat capacity is not 
a consequence of isomerizations. Therefore the general, 
two-peak shape of the heat capacity is not related to iso- 
merization, but is intrinsic to the potential itself and the 
presence of two types of particles organised in shells. 

The variation of the melting temperature with cluster 
size is shown in Fig. |S1 for the three {n,4}, {n,5} and 
{n,6} binary series, with 2 < n < 20. Here the melt- 
ing temperature was defined according to the Lindemann 
criterion, as the temperature where d sharply increases 
above 0.1. The variation of the computed melting point 
is qualitatively similar to that reported by Drocco and 



coworkers [13. In particular, the high values for {4,4}, 
{5,5}, {6,5} and {7,6} are reproduced. However, our re- 
sults quantitatively differ from those of Drocco et al. by 
a factor of at least 10. Our data, obtained with exactly 
the same Hamiltonian as in Ref. IT3, are in agreement 
with the computed melting points for the homogeneous 
clusters, which are comparable with previous results for 
the same systems 0. Even though the present temper- 
atures result from time averages in our MD simulations, 
their range matches the canonical temperature of our MC 
runs. Hence we believe that the data of Fig. |H|are reli- 
able, and that the differences from the results obtained 
by Drocco and coworkers are the consequence of an un- 
reported alternative choice of A in the Hamiltonian used 
in Ref. [H. 

Lastly, we discuss the behaviour of larger clusters and 
how melting evolves toward bulk. In Fig. [3 a) we show 
the heat capacities calculated from MC simulations for 
the cluster sizes n = 100, 200, and 500. The general 
trend is that the main peak near T ~ 0.04 becomes nar- 
rower and sharper, suggesting that melting is a first-order 
process rounded by finite size effects [s^l- However, the 
main heat capacity peak remains quite broad, so a more 
detailed study should probably be performed on even 
larger sizes, including a proper finite-size scaling anal- 
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FIG. 5; Disconnectivity graphs for the binary clusters of Fig. 0| The global minimum energy is shifted to zero for each size. 
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FIG. 6: Size evolution of the melting point of the {n,4}, 
{n, 5}, and {n, 6} series of binary Coulomb clusters as ob- 
tained from the Lindemann criterion. 



ysis. In addition, extra features in the heat capacity are 
present at low temperatures, even for the largest 500- 
particle cluster. Regular quenching confirms that these 
premelting effects are due to isomerizations involving dif- 
ferent arrangements of the inner shells. At large sizes, 
even though the global shape of the cluster is roughly 
circular, the radial ordering into shells is replaced by a 
trian gula r close-packed lattice, possibly with several de- 
fects [l^- Rotation of the various shells then becomes 
hindered, and the dynamics is reduced to more usual con- 
secutive isomerizations. 

Fig. EJa) supports the idea that a single, first-order 
phase transition should occur in the bulk limit, if it could 
be defined. Since that there is no unified theory of 2D 
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FIG. 7: (a) Canonical heat capacities (in units of the Boltz- 
mann constant) of some large homogeneous clusters; (b) melt- 
ing point of homogeneous clusters versus their radii. 



melting jSJ], we should probably be cautious in extrap- 
olating the heat capacities computed from Monte Carlo 
simulations. If we choose the peak of the heat capacity as 
the signature of melting, the melting point thus defined 
shows very pronounced variations with the cluster radius 
n^/^, as illustrated in fig. Cfb). 

The data for all sizes in the range 5 < n < 50 and 
for larger clusters (n = 60, 70, 80, 90, 100, 120, 140, 160, 
180, 200, 250, 300, 350, 400, 450, 500) are included in this 
graph. The strong size dependence for small systems is 
similar to that observed for the binary systems. However, 
in the range 50-500 the variation of the melting point 
seems less drastic than in the smaller size range, even 
though numerous sizes are obviously missing from the 
reported data. The bulk limit for the melting point is 
obtained to be 0.039. 

This observation should not obscure the possibility 
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that even large sizes can exhibit a multiple stage melt- 
ing. For example, significant premelting peaks have been 
reported in large atomic clusters The results ob- 

tained for n = 500 show a similar feature. Even though 
some premelting peaks are present in the caloric curves 
of large clusters, a major peak seems to dominate above 
about 100 ions. This result may suggest that the tran- 
sition between the finite size regime and bulk j35| takes 
place near this size. 



IV. DISCUSSION AND CONCLUSION 

In this work, we have investigated some dynamic and 
thermodynamic properties of finite 2D charged particles 
confined by a parabolic trap. Our extensive data on 
homogeneous and heterogeneous clusters were obtained 
from MD and MC simulations. These data were inter- 
preted by visualising the energy landscape through dis- 
connectivity graphs. Previous experimental and numeri- 
cal findings suggest that there are three types of motion 
for a particle: it can oscillate, rotate in a relatively regu- 
lar orbit, or make jumps between shells, possibly causing 
isomerizations. These various phenomena were also ob- 
served in the present work, but we found evidence that 
additional 'premelting' events due to structural transi- 
tions may also occur. The various stages of melting cre- 
ate ambiguities in defining the melting point, especially 
in small clusters. 

Firstly, rotational motion is facile when very small en- 
ergy barriers separate the global energy minimum from 
itself, as is apparent from the disconnectivity graphs 
for the 17- and 21-particle homogeneous clusters, or the 
{4,5} binary cluster. Internal rotation may give rise to 
significant jumps in the RMS bond length fiuctuation, 
sometimes above 0.1. Radial melting is then seen as 
yet another increase in 5 or other dynamical indicators 
at higher energy. Therefore, the Lindcmann criterion 
((5 ^ 0.1) is not fully reliable for detecting melting. 

Secondly, rotational melting has no thermodynamic 
consequence since only a single isomer is involved. The 
thermodynamics is sensitive to the presence of different 
isomers, and also to the intrinsic shape of the potential 
felt by the particles. Irrespective of the isomers, the heat 
capacity has a natural drop over a rather broad tempera- 
ture range, which marks the continuous transition toward 
the gas phase. Hence, even a single cluster isomer can 
have heat capacities with complex temperature depen- 
dence, as exemplified by the small {4,5} binary clusters. 

Thirdly, several isomers can give rise to preliminary 
peaks in the heat capacity, which may sometimes ob- 
scure the features associated with full radial melting. 
We found such a case for the 31-particle cluster inter- 
acting through dipolar forces, but very large clusters can 
also exhibit particularly strong premelting peaks. The 
structural transition occurring in the 31-particle clus- 
ter has the same origins as the fcc-icosahedral transition 
in the 3D 38-atom Lennard- Jones cluster |^ |53|: the 



second morphology has a much larger entropy than the 
global minimum. This phenomenon has important con- 
sequences for the constant energy dynamical behaviour, 
as this transition is not seen until sufficient energy to sur- 
mount the isomcrization barrier is present. Thus there 
is a finite energy range where MD simulations starting 
from the ground state structure are non-ergodic, which 
is reflected by a strong shift in the estimated melting 
temperature. 

Most of our results illustrate how a simple and unique 
picture of melting is hampered by finite size effects. 
Both the dynamical parameters and the thermodynami- 
cal curves contribute to explaining related but different 
aspects of melting in these 2D clusters. In particular, our 
work does not support a single parameter as a universal 
measure for the melting point that could be used unam- 
biguously for all clusters. Nevertheless, the landscape 
approach ^3 that has been followed here was found to 
provide useful insight for interpreting and understanding 
the various numerical data. 

The data gathered here for a broad range of sizes 
would seem to prevent a general classification on the 
phenomenology of melting to be made. For exam- 
ple, the choice of molecular dynamics or Monte Carlo 
method, and the microcanonical or canonical ensemble, 
can change the observed behavior significantly, particu- 
larly at small sizes where the differences between the two 
statistical ensembles are the largest. For instance, some 
structural transitions might not appear with conventional 
MD methods, due to the disconnected character of phase 
space at low total energies. 

However, three main categories of behavior can be in- 
ferred from our results. Especially in small clusters, a 
two-step orientational-then-radial melting has been ob- 
served, with clear dynamical signatures, but only a weak 
thcrmodynamical signature. Clusters exhibiting struc- 
tural transitions, on the other hand, may not involve any 
rotational melting but display a strong change in their 
caloric curve. Finally, a more conventional single-step 
melting is sometimes seen, as expected in larger clusters, 
and also for smaller 'magic' sizes. 

Our calculations indicate that significant finite size ef- 
fects remain in the thermodynamics of the largest clusters 
considered. However, above a few hundred ions, melting 
seems to be essentially a one-step process. This result 
is consistent with the expectation that the dynamics be- 
comes more and more dominated by localised processes, 
due to the less favourable rotational motion. This obser- 
vation seems to justify the use of the Lindcmann index 
or other dynamical parameters for larger systems. Ad- 
ditionally, the main heat capacity peak becomes larger, 
and its maximum value coincides with the onset of dis- 
order in the dynamical indicators. Even though premelt- 
ing features may remain noticeable, we expect the onset 
of radial melting to become closer and closer to orienta- 
tional melting, and eventually similar to a bulk first-order 
melting process. 
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